import random
import math
from IPython.core.display import SVG
import pyomo.environ as pyo
from pysat.solvers import Solver
from pysat.formula import CNF 
import py_svg_combinatorics as psc
from ipywidgets import widgets, HBox
from collections import Counter
from pprint import pprint
from random import randint
import numpy as np
from IPython.display import IFrame
import IPython
from copy import copy
import os
from pathlib import Path
np.set_printoptions(precision=2)
np.set_printoptions(suppress=True)
nbname = ''
try:
    nbname = __vsc_ipynb_file__
except:
    if 'COCALC_JUPYTER_FILENAME' in os.environ:
        nbname = os.environ['COCALC_JUPYTER_FILENAME']    
title_ = Path(nbname).stem.replace('-', '_').title()    
IFrame(f'https://discopal.ispras.ru/index.php?title=Hardprob/{title_}&useskin=cleanmonobook', width=1280, height=300)

def print_solution(m):
    for v in m.component_data_objects(pyo.Var):
        if v.value and v.value > 0:
            print(str(v), v.value)

m = n = np.random.randint(2, 3)

def generate_random_data_for_maxqp(n):
    m = n 
    # Линейные ограничения для допустимой области
    A = np.random.rand(m, n)
    b = 10*np.ones(m) #np.floor(1/np.random.rand(m))
    M = np.random.rand(m, n)
    Q = np.dot(M, M.T) # симметричная положительно-полуопределенная
    c = np.random.rand(n) 
    return A, Q, c, b


A, Q, c, b = generate_random_data_for_maxqp(n)
A, Q, c, b

def get_model_bbl(A, b, Q, c):
    m = pyo.ConcreteModel()

    m.m, m.n = A.shape

    # на всякий случай, возьмем с собой
    m.A = A
    m.Q = Q
    m.b = b
    m.c = c

    m.I = range(m.m)
#    m.J = range(m.n)

    
    upper_bound_x = math.ceil(max([b[i] / min([m.A[i, j] for j in m.I]) for i in m.I]))
    
    m.u = upper_bound_x #16 #256
    bin_len = math.ceil(math.log(m.u)/math.log(2))

    m.K = range(bin_len)

    m.x = pyo.Var(m.I, domain=pyo.NonNegativeIntegers)
    @m.Constraint(m.I)
    def не_больше(m, i):
        return sum(A[i, j] * m.x[j] for j in m.I) <= b[i]

    m.t_ik = pyo.Var(m.I, m.K, domain=pyo.Binary)

    @m.Constraint(m.I)
    def t_будет_бинарным_разложением_x(m, i):
        return m.x[i] == sum(2**k * m.t_ik[i, k] for k in m.K)

    # Бинарное произведение y_ikjl = t_ik * t_jl
    m.y_ikjl = pyo.Var(m.I, m.K, m.I, m.K, domain=pyo.Binary)

    @m.Constraint(m.I, m.K, m.I, m.K, [7, 8, 9, 10, 11, 13])
    def y_ikjl_constraint_A(m, i, k, j, l, c):
        if c == 7 and m.Q[i,j] < 0:
            return m.y_ikjl[i, k, j, l] <= m.t_ik[i, k]
        if c == 8 and m.Q[i,j] < 0:
            return m.y_ikjl[i, k, j, l] <= m.t_ik[j, l]
        if c == 9 and m.Q[i,j] > 0:
            return m.y_ikjl[i, k, j, l] >= m.t_ik[i, k] + m.t_ik[j, l] - 1
        if c == 10 and m.Q[i,j] > 0:
            return m.y_ikjl[i, k, j, l] >= 0
        if c == 11 and m.Q[i,j] != 0 and i < j:
            return m.y_ikjl[i, k, j, l] == m.y_ikjl[j, l, i, k]
        if c == 13 and m.Q[i,i] != 0 and k < l:
            return m.y_ikjl[i, k, i, l] == m.y_ikjl[i, l, i, k]
        return pyo.Constraint.Skip

    @m.Constraint(m.I, m.K, [12])
    def y_ikjl_constraint_B(m, i, k, c):
        if c == 12 and m.Q[i,i] != 0:
            return m.y_ikjl[i, k, i, k] == m.t_ik[i, k]
        return pyo.Constraint.Skip

    @m.Objective(sense=pyo.maximize)
    def сложная_побитовая_сумма(m):
        return sum( Q[i, j] * sum(2**(k+l) * m.y_ikjl[i, k, j, l] for k in m.K for l in m.K ) for i in m.I for j in m.I)

    return m

m = get_model_bbl(A, b, Q, c)
m.u
#m.pprint()

solver = pyo.SolverFactory('cbc')
solver.solve(m).write()

#print_solution(m)
print(m.x.extract_values())

[m.x[i].value for i in m.I]

import picos as pic
import matplotlib.pyplot as plt
import scipy
import cvxopt as cvx
import cvxpy as cp

x = cp.Variable()
y = cp.Variable()

constraints =  [x + y == 1,
                x - y >= 1]

obj = cp.Minimize((x - y)**2)

prob = cp.Problem(obj, constraints)
prob.solve() 
prob.status, prob.value, x.value, y.value

Q = np.array([[0, 1], [1, 0]])
Q

y = cp.Variable(2, integer=True)
obj = cp.quad_form(y, Q)
constraits = [y >= -2, 
              y <=  2]
prob = cp.Problem(cp.Maximize(obj), constraits)
# prob = cp.Problem(cp.Minimize(obj), constraits)
print(prob)
prob.solve() #solver=cp.SCIP

n_points = 5
u = np.linspace(-2, 2, n_points)
u

x1, x2 = np.meshgrid(u, u)
x1, x2

all_x = list(np.array([x1,x2]) for x1, x2 in zip(np.ravel(x1), np.ravel(x2)))
all_x

Z =  np.array([ np.dot(x.T, np.dot(Q, x)) for x in all_x])
Z = Z.reshape(x1.shape)
Z

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection='3d')
plt.figure()
ax.plot_surface(x1, x2, Z)

A, _, c, b = generate_random_data_for_maxqp(2)
A, Q, c, b

m = get_model_bbl(A, b, Q, c)
solver.solve(m).write()

[m.x[j].value for j in m.x]